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TEST PARTICLE PROPAGATION IN 


MAGNETOSTATIC TURBULENCE 
I. FAILURE OF THE DIFFUSION APPROXIMATION 

Abstract 

Test particle propagation in magnetostatic turbulence with a strong mean 
field component is considered. The equation which governs the quasi-linear 
approximation to the ensemble and gyro-phase averaged one-body probability 
distribution function is constructed from first principles. This derived equation 
(the quasi-linear diabatic equation) is subjected to a thorough investigation in 
order to calculate the possible limitations of the quasi-linear approximation. 

It is shown that the reduction of this equation to a standard diffusion equation in 
the Markovian limit can be accomplished through the application of the ’’adiabatic” 
approximation. It has been shown that this standard diffusion equation is identi- 
cal to that obtained through the assumption of a Fokker-Planck equation for the 
probability distribution function. In the presence of the strong mean magnetic 
field, the reduction to the Markovian limit is shown to be invalid. Numerical 
solutions of the (integrodifferential) quasi-linear diabatic equation are obtained 
using a simple axisymmetric slab model of the turbulent field for (i) narrow 
parallel beam injection, (ii) broad parallel beam injection, and (iii) narrow cross- 
field beam injection. A numerical solution of the standard diffusion equation in 
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the Markovian limit is obtained for the narrow parallel beam injection. Com- 
parison of the diabatic and adiabatic results explicitly demonstrates the failure 
of the Markovian description of the probability distribution function. This failure 
is discussed in terms of an appropriate mode ("Laplace-mode") expansion of 
these solutions. For parallel beam Injection, the relaxation to isotropy is shown 
to proceed slowly according to ln(time)/time, in contrast to the familiar expo- 
nential relaxation usually associated with diffusive behavior. The relaxation of 
a parallel beam to isotropy in the Markovian approximation is shown to proceed 
with a ziero exponential decay rate; i.e., isotropy is never reached. Through the 
use of a linear time-scale extension the failure of the adiabatic approximation, 
which leads to the Markovian limit, is shown to be due to mixing of the relaxation 
and interaction time scales in the presence of the strong mean field. 



I. INTRODUCTION 


This is the first in a series of papers in which the propagation of a charged 
test particle in magnetostatic turbulence is considered. The central issue in this 
work is the presence of a mean magnetic field which is strong enough to affect 
the collision process between the charged particle and the turbulent field. For 
the purpose of this study, the magnetic field is considered ’’strong” if the Larmor 
radius of the particle in the mean field is comparable to the two-point correlation 
length associated with the turbulent field. If electron- ion collisions in a plasma 
were being considered here, the field would be considered ’’strong” if the Larmor 
radius were comparable to the Debye length. In this situation some standard 
methods of kinetic theory which are suitable in studies of neutral gases or non- 
magnetized plasmas, become inapplicable. The failure of these methods will be 
demonstrated here, and alternatives will be presented. 

Test particle propagation in magnetostatic turbulence can be considered an 
idealization of the behavior of a particle in a hot, turbulent MHD plasma in which 
the effects of particle-particle collisions are negligible compared to the effects 
of particle-wave collisions. If the speed of the particle is high compared to the 
MHD wave propagation speeds, then the motion of the particle is dominated by 
the magnetic field, and, as interesting as the effects of the electric field are, 
they are nevertheless, secondary. The full, self consistent, treatment of the 
plasma problem is not attempted here; if it were, then the work being presented 
here would be a necessary part of that attempt. Alternatively, sufficient 
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measurements of the field can be relied upon to test this portion of the complete 
treatment. The solar wind plasma which fills interplanetary space is an excel- 
lent candidate for this type of study; it is a dilute, high temperature, turbulent 
MHD plasma whose particle and field properties have been measured extensively.^ 

The interplanetary plasma has been observed to have an extremely long, 
power law tail in its energy distribution. The particles which make up this tail 
are called cosmic rays. They have their source predominantly at the Sun in the 
several tens of Mev/particle energy range, outside the solar system for higher 
energies, and perhaps, even outside the Milky Way galaxy at the highest observed 
energies 10 eV/particle).^’ ^ The speeds of these particles are usually high 
compared to the ts^ical MHD wave speeds found in the interplanetary or inter- 
stellar plasma. Thus, the assumption of magnetostatic turbulence is a very 
good approximation for these particles. In addition, the energy density of these 
particles is so low, they almost always have negligible affect on the field through 
which they move in any physical system which has dimension less than that of an 
entire galaxy like the Milky Way. ^ As might be expected, there is a long history 
of the treatment of the cosmic ray propagation problem in the magnetostatic, 
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test particle limit. 

The earliest attempts to describe the spatial transport of cosmic rays were 
invariably based on diffusion equations for the cosmic ray density with spatial 
and energy convection terms added when it was deemed necessary. ® With the 
advent of space exploration, and in particular, with the advent of detailed 
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measurements of the interplanetary magnetic field, it became possible to investi- 
gate the theoretical foundations of these assumed diffusion equations. In the first 
efforts in this direction, kinetic theories which were known to be applicable in 
other situations were borrowed and assumed to apply to the cosmic ray prob- 
ability distribution function. Klimas,^ and independently, Gleeson and Axford^ 
adopted the Boltzman equation and through moment expansions of the distribution 
function, constructed transport equations which were similar to those in use pre- 
viously except that the transport coefficients all became interrelated through an 
effective mean free path for scattering on the magnetic inhomogenieties. This 
effective mean free path came from the Boltzman collision integral which was 
assumed to represent the wave-particle interaction. Application of these theories 
was accomplished through phenomenological adjustments of the effective mean 
free path. The relationship between this mean free path and the actual inter- 
action mechanism remained vague. 

Jokippi,® and independently Sturrock,^ adopted the Fokker-Planck equation 
to describe the phase-space propagation of the cosmic ray probability distribu- 
tion function. These theories were also convection-diffusion theories in phase- 
space, and through moment expansions, could be used to construct convection- 
diffusion transport theories in configuration space which were again similar to 
those in use previously. The major significance of these theories was that the 
transport coefficients which appeared in them could be directly related to meas- 
urable properties of the interplanetary field; specifically, the mean field and the 
two-point correlation tensor associated with the turbulent field. 
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These results were reassuring in a certain sense; no matter what kinetic 


basis was considered, a convection-diffusion transport theory seemed to result. 
However, what had really been accomplished was to make the assumption of a 
diffusion process somewhat less visible. Both the Boltzman, and the Fokker- 
Planck equations can be derived from first principles under appropriate condi- 
tions and in their respective realms of applications.^^The adiabatic approximation, 
or its equivalent, is a necessary part of these derivations. But, the imposition of 
the adiabatic approximation is equivalent to the assumption of a Markov chain for 
the relevant collision process. If the adiabatic approximation can be applied, then 
the particle motion is well approximated by a random walk process, and the macro- 
scopic transport description of the fluid in question is necessarily governed by a 
diffusion equation with convective phenomena possibly included. 

The first attempt to derive the kinetic theory of test particle propagation in 
magnetostatic turbulence from first principles was made by Hall and Sturrock.'^^ 
They managed to regain the Fokker-Planck equations which had been assumed by 
Jokipii and Sturrock, but, in doing so, they had to apply two approximations; first 
the quasi-linear approximation, and then, in fact, the adiabatic approximation. 

In the quasi-linear approximation, the impulse imparted to a particle during an 
interaction with the turbulent field is assumed small compared to the momentum 
of the particle. Thus, this impulse is calculated on the basis of an undeviated 
trajectory through the interaction. Implicit in this picture, is the assumption 
that the interaction between the particle and the turbulent field can be characterized 
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as a weak "collision” of short duration so that the net impulse imparted to the 
particle remains small. In keeping with this assumed weak coupling of short 
duration, it seems reasonable to assume that the distribution function should 
evolve on a time scale which is much larger than the duration of a collision. 

Thus, it seems reasonable to make the adiabatic approximation, in which the 
evolution of the distribution function during a collision is ignored. This picture 
is correct if the mean magnetic field is not strong, or alternatively, if the par- 
ticle energy is high.^^ However, Klimas and Sandri,^^ using a special isotropic 
model of the magnetostatic turbulence and a spherical harmonic expansion of the 
distribution function, showed that this picture becomes incorrect when the field 
becomes strong. Klimas and Sandri showed that, within their model, when the 
mean magnetic field is strong, then it is inconsistent to apply both the quasi- 
linear and adiabatic approximations. This inconsistency arises because in the 
presence of the strong mean magnetic field, the undeviated, or zero'th order 
trajectory, is a helix and therefore progress of the particle in space in the direc- 
tion of the mean magnetic field is governed by the parallel component of its 
velocity which is a constant of the particle motion. If the parallel velocity is 
made arbitrarily small, then the particle becomes quasi-trapped in its inter- 
action with the random component of the field, and one, or both, of the assump- 
tions on the smallness of the impulse imparted to the particle during the inter- 
action, as well as on the clear separation between the interaction and evolution 
times scales, breaks down. 
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Recently, Golstein, Klimas and Sandri^** have shown that, within the quasi- 
linear adiabatic system of approximations (or equivalently, in the Fokker-Planck 
equation), when the parallel component of the particle velocity is zero, then with 
few exceptions, the calculated strength of the interaction with the random field 
becomes infinite. They have also shown that this divergence can be directly 
attributed to the physical phenomenon of mirroring. Clearly, if the impulse im- 
parted to a particle through its interaction with the random field results in mir- 
roring of the particle, then this impulse in the parallel direction cannot be con- 
sidered small. In the approximation scheme being discussed here, in which the 
particle-wave interaction is considered asymptotically small, the finite impulse 
due to mirroring appears as a divergence in the theor}'. 

One of the few exceptions to the appearance of a divergence in the wave- 
particle interaction strength as calculated in the quasi-linear adiabatic approx- 
imation, is found in the slab model of the random magnetic field which we con- 
sider in this series of papers. In this model of the field, the interaction strength 
is calculated to be zero at zero parallel velocity. Because first order mirror- 
ing is impossible in this model of the field, the divergence due to mirroring 
vanishes. However, the interaction and evolution time scales still remain mixed, 
and a more subtle failure of the quasi-linear adiabatic approximation scheme 
will be explicitly demonstrated with the use of a linear time-scale extension in 
Section IV. It is the general failure of the quasi-linear adiabatic approximation 
due to the mixing of the interaction and evolution time scales which is of concern 
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to US here, more than any specific field-model dependent manifestiation of this 
failure. 

Various attempts have been made to construct non-linear theories of test 
particle propagation in which the undeviated particle trajectory is replaced by 
one which contains the effects of the wave-particle interaction being calculated,^^ 
In these theories, it is assumed that the particle propagation in the random field 
can be described through a diffusion process, and the undeviated trajectory is 
replaced by a diffusing trajectory. The amount of diffusion in the trajectory is 
computed so that the diffusion coefficient which is calculated using that trajectory 
is self consistently calculated. In this case, it becomes impossible for a particle 
trajectory to remain trapped, and with this modification of the quasi-linear ap- 
proximation, the adiabatic approximation can be formally retained. The final 
result, then, is a Markovian, or diffusive, descriptio’’. of the particle transport 
which is consistent with the diffusive modification of the undeviated trajectory. 

In this method, the adiabatic approximation is retained, and the quasi-linear 
approximation is modified to allow that choice. The reasons for this choice are 
unclear. In particular, we believe that the predictions of the quasi-linear theory, 
with no further approximations, have never been calculated. It has not been clear 
what, if anything, is incorrect in the quasi-linear theory. If the quasi-linear 
theory does fail , then modifications of that theory should be made on the basis 
of that failure. With this approach, perhaps a kinetic theory can be constructed 
which is free from initial prejudices on the ultimate result of that construction. 
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In this series of papers, we determine the predictions of the quasi-linear 
theory in a very simple model of the turbulent magnetic field which will be fully 
described later. We find that the adiabatic approximation to the quasi-linear 
theory does not make sense. We further find an alternative kinetic approxima- 
tion which works very well, artd is not governed by a diffusion equation. We 
conclude that the quasi-linear theory does fail as a leading approximation in a 
systematic expansion scheme, and we demonstrate why this is so. 

In this first paper, we construct an equation from first principles for the 
quasi-linear approximation to the ensemble and gyrophase averaged probability 
distribution function. This equation is a velocity space diffusion equation which 
is non-local in time; it is integrodifferential in time. We demonstrate that appli- 
cation of the adiabatic approximation to this equation leads to a Markovian de- 
scription of the probability distribution function which is governed by an ordinary 
velocity space diffusion equation. This diffusion equation can also be obtained 
through the assumption of a Fokker-Planck equation for the probability distribu- 
tion function.^^ Through both analytic and numerical considerations, we demon- 
strate the failure of the Markovian, or diffusive, description of the quasi-linear 
probability distribution function. Thus, we conclude that the presence of a strong 
mean magnetic field can preclude the application of a Fokker-Planck or diffusion 
equation to the propagation problem for charged particles in plasma turbulence. 

We proceed with a study of the properties of the quasi-linear solutions with 
no further approximations in this paper. Numerical solutions of the quasi-linear. 
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and quasi-linear adiabatic equations are obtained, and compared, for an axi- 
symmetric, slab model of the plasma turbulence with an exponential two-point 
correlation function. In this model, the random component of the field is orthog- 
onal to the mean field direction, it is stationary in time, and a function of the 
spatial coordinate in the direction of the mean field only. In addition, the two- 
point correlation tensor associated with the random component of the field is 
assumed cylindrically symmetric about the mean field direction. With this model 
the quasi-linear adiabatic theory predicts that the parallel velocity of a particle 
can never reverse itself as a result of the interaction with the turbulent field. 

This prediction is in strong disagreement with the numerical solutions of the 
quasi-linear equation without the adiabatic approximation (the quasi-linear di- 
abatic equation). 

We conclude this paper with a discussion of a mode expansion of the solutions 
of the quasi-linear diabatic equation. We introduce a "Laplace-mode" analysis 
which is analogous to the standard modal analysis often used in plasma physics. 
We show that the long time evolution of the quasi-linear solutions can be under- 
stood in terms of these Laplace-modes. In addition, we find that the reduction 
of these solutions to their adiabatic limit is easily expressed with the Laplace- 
modes, and that the failure of this limit becomes apparent within this picture. 

For example, the injection of a beam of particles in the direction of the mean 
magnetic field is considered in some detail. As pointed out above, in the quasi- 
linear adiabatic approximation, no particles ever reverse their directions, and 



so, the beam never relaxes to isotropy. In the quasi-linear diabatic approxima- 
tion, the beam does relax to isotropy, but very slowly. From the Laplace-mode 
analysis we find that in this case the relaxation to isotropy goes like ln(time)/time, 
in contrast to the exponential relaxation usually associated with diffusive behavior. 
(In the adiabatic case, the diffusive relaxation (tf the beam is exponential, however, 
with a zero decay rate. Thus, the beam never relaxes.) In the adiabatic limit, 
we find that the Laplace-modes become doubly degenerate, and discontinuous func- 
tions of the parallel component of the particle velocity. The failure of the adiabatic 
limit follows from these discontinuities in the Laplace-modes. In the next paper 
in this series we introduce a new kinetic approximation which is constructed spe- 
cifically to remove these discontinuities in the Laplace-modes. 

n. THE BASIC EQUATIONS 

We imagine an ensemble of stochastic, stationary magnetic fields. A mean 
field, <(|^, which is an ensemble average is assumed, and the field in any ensemble 
representative is represented by B = *(b)> + B’. The "random field," B’, is 
assumed to obey <^B’)> = 0, or equivalently, it is assumed that = <^B^. 

In the absence of two-body or higher order particle correlations the one-body 
probability distribution function, F(x, p , t), in any one of the ensemble represent- 
atives obeys the Liouville equation, 

|^ + KF + £F + t 7C'F= 0 I.i 

3r 
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The dimensionless version of the Liouville equation given here has been obtained 
by setting r = Wq t , where is the Larmor frequency in the mean field strength, 
and by measuring lengths in units of the gyro-radius, r^, in the mean field 
strength. The parameter, t] , is a measure of the strength of the random field 
compared to the strength of the mean field; it is defined by. 


B' 

_ rms 

V =■ 


<B> 


1.2 


and will be assumed small. The linear, first order partial differential operators, 
K, C , and, , are given by. 


K = p • 




1.3 


and. 


£=(PkP)-^=-P-S2-^ 

Bp Bp 


£' = (pxp') ■ _ p .Q' . _1 


1.4 


1.5 


Bp ^P 

A ^ 

where P is a unit vector in the direction of the mean field, p is a unit vector in 
the direction of the particle momentum, p' = B'/B' , and B/Bp is defined by 


B \ ^ . B 

= P(.hi -Pi Pj) 




1.6 


The skew-symmetric tensors, £2 and ft' , are defined by 

j ” ^ijk^k • I^^ter, the integral or correlation length associated with the two- 
point correlation function for the random field will be introduced. In this work 
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this length will be assumed comparable to the particle gyro-radius, and as a 
matter of convenience we set 



1.7 


e 

We also introduce an ensemble averaged probability distribution function, 
f(p ,T) which is assumed independent of position, x , and the random prob- 

ability distribution function through F(x,p,r) = f(p, r) + F'(x,p,r). We further 
assume <(F'^= 0, or <^f)>= f. Klimas and Sandri,^® using a technique developed 
by Kaufman,^’ have shown that, if F'(r = 0) = 0, then 


+ £f = [1 - 7] <£'G>]"' <£'G£'> f 1.8 

where, G = [B/^r + Ji] \ is the Green's integral operator, with the total Hamil- 
tonian operator, 34 = K + £ + ?]£', for the generator of the particle motion. The 
Green's integral operator can be written, 

G= [1 +7]Go£']'^ Gg 1.9 

where, G^ = [B/Br+ 34 ^]"^ , with the zero'th order Hamiltonian operator, 34 ^ = £ + k, 
for the generator of the particle motion. (We assume that the mean field is uni- 
form in space, and then, 34 q generates the well known helical particle trajectories 
in a uniform magnetic field.) With the use of equation 1.9, we expand equation 1.8 
in powers of t ] Gq£' to obtain, 

+ £ f = 7,2 <C G„C'> f t 0(V) 1.10 
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If this expansion is valid, then we say that we are in the ’’weak coupling” regime, 


or that we have made the weak coupling assumption. A truncation of this expansion 
beyond the second order term which is exhibited explicitly on the right hand side 
gives an equation for the quasi-linear approximation to the ensemble averaged 
probability distribution function. A major purpose of this series of papers is to 
present a class of situations in which the formal order of the terms which have 
been dropped in this truncation can be shown to be incorrect when the quasi-linear 
approximation to the probability distribution function is assumed a valid approxi- 
mation and is used to evaluate the actual order of the correction terms in this 
approximation scheme. 

III. THE QUASI-LINEAR TRUNCATION IN THE AXISYMMETRIC 

SLAB FIELD MODEL 

We adopt the quasi-linear truncation of equation 1. 10, and develop the specific 
form that that equation takes on in an "axi symmetric slab” model of the random 
magnetic field. In this model, the random field is orthogonal to the mean field, 
and a function only of the spatial coordinate which lies along the direction of the 
mean field. In addition, the random field takes on any direction in a plane orthog- 
onal to the mean field with equal probability from ensemble representative to 
representative. 

The particle energy is a constant of the motion in the static magnetic field. 
Since we have already assumed that the ensemble averaged probability distribu- 
tion function is independent of position, we see that the phase space relevant to 
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our problem is the surface of a sphere at constant energy in the particle momen- 
tum space. Actually, through our dimensional analysis which lead to the dimen- 
sionless Liouville equation, we have reduced our problem to that of studying the 
motion of the particles on the unit sphere. Thus, it is most convenient to intro- 
duce a polar coordinate system, {6 ,<p) in the particle momentum space, in which 
the polar coordinate, 8 , is the pitch angle of the particle relative to the mean 
field, and 0 is the gyro-phase angle of the particle measured about the mean 
field, hi addition, we introduce, /x = cos 6. The differential operators, I and £' , 
are given in terms of these variables by, 



1.11 


= (p .Q.p') ^ _iLJl (p-p') ^ 


V 1 _ ^2 


30 




1.12 


The gyro-phase average of any quantity, Q ( 0 ), is defined by, 


Q = — f d0 Q(0) 
*'0 


1.13 


Through the substitution of equations 1. 11 and 1.12 into the truncated equation 1. 10, 
and further, through gyro-phase averaging the entire equation, we find. 
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~9f (m. t) _ 


= - 77^ 




d\ I- p • Q (O P' (.^)>- a • P 

OfJi OfJ. 

1,14 


772^1* d0 f d\ A p • n (z) e p . ^ . \ ^ 

Jq Jq ^P- \1 - ij.y 


where we have adopted the symbol, z, for the spatial coordinate along the mean 
field, and where the zero'th order streaming operator, exp(-J^P \) acts on any 


function, A{z,p), as follows: 


with, 
in which. 


_'iy \ 

e ° A(z, p) = A(z - yu\) C(\) • p) 


C (X.) = P+ IN cos X-Q sinX. 


P. . = /3. 
N. . = S.. 

XJ X3 


1.15 

1.16 

1.17 

1.18 


and Q is as defined previously. We assume that the magnetostatic turbulence is 
homogeneous and introduce the two-point correlation tensor and its axisymmetric 
slab model through, 

</3!(z) + 0> " Rij(0 = N,.R(0. 1.19 

with the correlation function, R(^ ), an even function of 1. Through the substitu- 
tion of equations 1.15 thru 1.18 into equation 1.14, we find the second term on the 
right side of that equation is identically zero, and the remainder of that equation 
reduces to. 
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(1 . ^2) 

dr d/i 


r 

•>0 


dX. K(jj., X) 


3f Qu, r - X) 

'dfJL 


1.20 


where the kernel function, 3{(Ai , X), is given by, 

X) = R(/iX) cos X. 1.21 

Thus, in the axisymmetric turbulence, the gyrotropic part of the probability dis- 
tribution function evolves independently. 

Equation 1,20 can be written in the renewal form, 


f (jXy T) = f O, 0) + 77 ^ “ (1 - 

Ofl 


■X 


f dX f ds KQi, s ) iL (^^'^-^ h.22 


and, with the introduction of 


N(Ai, r) 


i M 

d/x' 

1 


fC^'. r) 


1.23 


the equivalent equations, 


(1 - ^2) f dk K(js. K) -r - 

Jo -bp? 


1.24 


and, 


N(/x, r) = NO, 0) + 772^1 - ;x2) f dX f ds 3((/x* s) ^ 1.25 

Jq Jq B/j.2 


can be constructed. All of these equations will be used in the following to under- 
stand the behavior of f (/x , r). 
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IV. THE MARKOVIAN LIMIT 


Equation 1.20 can be characterized as a non-local diffusion equation; it is 
integrodifferential in time. This equation can be further reduced to a standard 
diffusion equation through the application of the "adiabatic" approximation. 

This diffusion equation is identical to the diffusion equation which follows from 
the assumption that the ensemble averaged probability distribution function obeys 
a Fokker-Planck equation. Since the Fokker-Planck equation has a Markovian 
process as a fundamental assumption, the adiabatic approximation reduces the 
non-local diffusion equation, which has an obvious memory of the past, and is 
non-Markovian, to its Markovian approximation. 

We will demonstrate in this section that a necessary requirement for the 
reduction of the non-local equation to the local one, is a clear cut separation 
between the presumed fast time scale during which a "collision" with the turbulent 
magnetic field takes place, and the much slower time scale over which the distri- 
bution function is assumed to evolve. In fact, this separation of the two time 
scales does not exist, and the Markovian approximation can not consistently be 
made. We will first present a brief intuitive derivation of the adiabatic approxi- 
mation and then a more rigorous derivation based on the time-scale extension 
method.^° We will see that, although the extension method is able to produce a 
imiform expansion of f (/i , r) in time, the expansion is still non-uniform in 
The non-uniformity in p- can be traced directly to the mixing of the interaction 
and evolution time scales in the m -domain in which the expansion is non-uniform. 
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(a) Preliminary Construction 


The intuitive argument which is often used to justify the adiabatic approxi- 
mation proceeds as follows: One can see from equation 1.20 that the time rate of 
change of f (Ai-.r) is small if 77 is small. Therefore replace f(/^, '7'-^) in the inte- 
grand of that equation by f . '^) when r is large and ^ is not. Rely on the pre- 
sumed short range of the kernel to prevent contributions to the integral for large 
i.e., assume that = 0, if A- >>1. Then, for '^>,> 1 , 

1.26 

where, 

^q(J^) ~ (1 - I d\ 3((At, k). j^27 

"0 

Equation 1.26 gives the Markovian, approximation to f(/^,T). This equation 
has been obtained by assuming that f (a^.t) does not evolve in time over the short 
period of time during which (/x A) - 0 with increasing X . From equation 1.21, 
the short range of the kernel should be provided by the short range of the corre- 
lation function in the random field, but the range of the correlation function in k 
can be made arbitrarily large as a*- is allowed to take on values arbitrarily close 
to zero. It is this mixing of the evolution and interaction time scales for those 
particle trajectories which are quasi-trapped in the interaction with the random 
field due to their small parallel velocity, a^, which leads us to suspect the Mar- 
kovian description off (a ,r). 
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(b) Linear Time- Scale Construction 

Through the use of a linear time scale extension^® we are able to obtain 
equation 1.26 as a leading result in a systematic expansion procedure which also 
yields a higher order correction to f (/i ,r). We will demonstrate here that this 
correction becomes unbounded for large times when /j- = 0, thus providing evi- 
dence for the failure of the adiabatic, or Markovian, approximation to 
The extension of equation 1.20 which we consider here is given by, 

^ -f T]2 3(/x, Tg, 1-2) = ^ (1 -/i.2) J d\K(p., k) Tq-K 


Br, 


The "restricted trajectory" is characterized by Tq = r and T2 = and on this 


restricted trajectory, we require 

Q(jjl, r, rj^r) = i (jm, t) 


1.29 


Thus, on the restricted trajectory, equation 1.28 reduces to equation 1.20. We 
further introduce, 

30. 'Tq. '^ 2 ) = 3qO. 7-g, r^) + 7]^ \(p., Tg, T 2 ) + 0{Tf) 1.30 

and expand equation 1.28 in powers of . By equating coefficients of powers of 
Tf we find, 

B3 q (jS, Tg , ) 


Br,, 


= 0 


1.31 


and 

B3j O’ ”^0 ’ "^2 ^ ^3q O’ '^0 ’ ”^2 ^ B 

Br„ Br„ 3^ 


+ = ^ (1 - IJ?) I ° dr 3<0a, 

r\ 


^3oO. '^0 '^2) 

B/x 


1.32 
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From equation 1.31, 


'^0* '^2) = ^o(^’ "^2) 


1.33 


and then, from equation 1.32, 




BSqO, 0, Tj) 




i m 

d>.J{0u, \) 


33q0u, 0, Tj) 

3/x 


^ (1-^2) 


f” , V, ,B3,(^,0,r,) 

J J ds s) X.34 


B/i 


We remove the secular growth of 3„ by setting 


aa„&., 0, r,) 3 f” 0, T,) 


Br„ 


ByU. 


I 


(1 - ^2) ci\ 3(0, 


B^ 


1.35 


and then, along the restricted trajectory, we find, 




Br 


= 7?2 


B/x 




B|L(. 


and, 


r) - i^(pL, 0) - — (1-^2) I d>w j ds K(ji, s) 


0 ''k 


t) 

B/x 


1.36 


1.37 


The behavior of Iq (/x, r), as r depends on the properties of (/x) = 

(1 - /U.2) (jj.), where IHq (/x) is the zeroth moment of the kernel. With the intro- 

duction of the power spectral density, 
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dis R(s) cos ws 


1.38 


we find that, 


P(^)^ 


r 

•'A 


ffloOu) = 




1.39 


We will assume the best possible case; i.e., we will assume that P(<^) is non- 
negative, and has no zeroe*i in the range 1 < oj < “ , but, of course P(°° ) = 0. In 
fact we must have , 



“* CO) 


1.40 


so that the total power in the random field remains finite. Thus, it is very gen- 
erally true that 

\Cc4)-0 1.41 


and, in the "best possible case" being considered here, HIq (/i.) has no other zeroes. 
Thus, J9p (fj,) is non-negative, and has zeroes at m = 0, ±1. 

The consequences of the zero in 19^ (p.) at fj. - 0 are developed thoroughly in 
Appendix B, but the discussion there depends on the normal mode expansion of 
the probability distribution function which will be constructed in section VI. In 
the following paragraph we give a brief description of the results of Appendix B 
and discuss the consequences of these results in the expansion being attempted 
here. 
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The solutions to equation 1,36 must approach a final steady state in which 
becomes independent of everywhere that (/^) ^ 0,, Because of the 
zeroes in at /x = ± 1 , equation 1.36 conserves the total probability contained 
in the domain, - 1 < /x < 1 . Because of the zero in (^) at ^tx = 0, we show in 
Appendix B that the total probabilities in the half-domains, -1 < /u. < 0, and 
0 </x < 1, are also individually conserved. Generally speaking, Io(/x,r) approaches 
a final state in which a discontinuity at /u. = 0 appears. Thus, 


BfpCO. -r) 

B/x 


~ 00 (7- -» 00) 


1.42 


From equation 1,37, we see therefore, that ,r) must also become unbounded 
at ^ = 0 with increasing time. The inability of equation 1.36 to allow the propaga- 
tion of probability through = 0, leads to a discontinuous f^{^t,r) at/x = 0, which 
in turn, leads to divergences in higher order corrections to Fo(/a,r) that invalidate 
the entire expansion procedure. 


The usual Markovian approximation to the non-local diffusion equation, equa- 
tion 1.20, cannot be made here. The solution of equation 1.20 is not adequately 
described by the ordinary diffusion equation, equation 1.36. The failure of the 
Markovian approximation can be traced to mixing of the time scale over which 
the probability distribution function evolves with the time scale during which an 
interaction with the random field occurs. This mixing is a fundamental problem 
in the quasi-linear expansion scheme in the presence of a "strong’' mean mag- 
netic field. In this expansion scheme the particles are carried through an 
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interaction with the random field along the undisturbed helical trajectories in the 
uniform mean field. Those particle trajectories which have small components of 
velocity along the mean field direction become effectively trapped in the inter- 
action with the random field. The interaction time, in these cases, becomes 
arbitrarily long in contrast with our assumption that it is short compared to the 
time over which evolves. Thus, we come to the important conclusion, 

that in the presence of a strong mean magnetic field, the use of the quasi-linear 
expansion scheme can preclude the Markovian, or diffusive, description of the 
probability distribution function. Of equal importance is the converse, that the 
Markovian, or diffusive, description of the probability distribution function can- 
not be used to judge the validity of the quasi-linear expansion scheme itself. 

V. NUMERICAL SOLUTIONS 

In this section we present the results of numerical integration of equation 
1.20 for several important initial conditions. The Markovian approximation to 
the quasi-linear solution, which was discussed in the previous section, has also 
been obtained numerically. A comparison of these two solutions, with the same 
initial conditions, will be given here. 

Numerical integration of equation 1.20 was found to be of limited use due to 
the large amount of computer time and core necessary to compute the convolution 
integral in that equation at each integration step. With the introduction of a spe- 
cial form for the correlation function. 
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1.43 


R(0 = e‘'^' 


we found that equation 1.24 could be considerably reduced to a system of three 
coupled partial differential equations (I.46-I.48) so that numerical integration 
became feasible. (We have investigated the consequences of using a double expo- 
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nential correlation function of the form suggested by Chernov in order to satisfy 
the requirement introduced by Khintchine that the correlation function have a 
zero first derivative at its origin, but have not found any qualitative differences.) 

With the exponential correlation function given by equation 1.43, equation 1.24 
is equivalent to. 


3N 2 
- 

Br 


1.44 


and. 


Bg I I .. 2\ , 

~+ l/^l g = (1 - - h 

3-^ dfj? 


Br 


+ \i^\ h = g 


1.45 


1.46 


where , 


and. 


gOit, r) ^ (1 f d\ e~ ^ cos \ 

Jo 'dfj? 




f 


'd/jr 


1.47 


1.48 


with the following boundary-initial conditions: 
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g(^, 0) = 0 h(/x, 0) = 0 

g(±l, r) = 0 h(±l, T) = 0 


1.49 


We treat ,0) as a given function of and notice from equation 1.24 that 
N(±l,r) is independent of time. From equation 1.23, we see that N(-1 ,t') = 0, and, 
since equation 1.24 is homogeneous, we choose N(+1 ,t) = 1 without loss of gen- 
erality. With the exponential correlation function, equation 1.24 is also equiva- 
lent to. 


Br 


B^N 1 BN 

+ ^ + (1 

Br2 Br 


-d/J 


= U1 (1 - 


B^N 

B/U.2 


1.50 


with, 


BN(/x, 0) _ n B^NCa, 0) _ ^2 _ ^, 2 ^ B^NQx, 0) 

Bt- 37-2 B/U.2 


1.51 


It is instructive to consider the structure of equation 1.50 in order to establish 
some contact with more standard equations of mathematical physics, and to sug- 
gest some qualitative features of its solutions. For this purpose, it is convenient 
to rewrite equation 1.50 as the equivalent pair. 


B2n b^n 


bn 


-V2(/x) + v(^)^+^^(^)N = C0u, r) 


'd/j? Br2 


Br 


I.52a 


BC . B2n 


B'^ B^2 


L52b 
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where, 


V^O) = 77^(1 - ^2) ^(/i) = 1 VO) vO) I.53a 

V O) = 2 l/^l (/^) = 1 + I-S3b 

If C{fM,r) is set to zero in equation I.52a, it reduces to a general dissipative Klein- 
Gordon wave equation.^ ^ The coupling between equations I.52a and I.52b, when 
C{ii,r) 4 0, exhibits an interaction between hyperbolic (wave-like) and parabolic 
(diffusive-like) behaviors in the single equation 1.50. There are four possible 
special cases of the homogeneous equation which depend on there being special 
relationships between the coefficients in equations 1.53: 

(a) The scalar wave equation 

V = X = 0; V = 

(b) The dissipative scalar wave equation 

X = 0; V = V - 0-/6 1.54 

(c) The Klein-Gordon equation 

V = c; V = 0 

(d) The telegraph equation 

V = 2(a + S); = 4a§; (a, S) real 

In case (d) the solution N(/u.,r), is analogous to the propagation of a voltage along 
a cable, where V = (LC)"^ , a = R, S = G/2C, R is the resistance, L is the self 
conductance, G is the capacitance, G Is the leakage conductance, o- is the con- 
ductivity, p is the permitivlty, and e is the dielectric. 
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Notice, if we were to neglect the first three terms of equation I.52a, so that 
C(/Li,r) = we would obtain, 

bn ^ (» - ^ 1.55 

(1 + p?) 2/J? 

which is just the diffusion equation which generates the adiabatic approximation 
to the probability distribution function in the special case being considered here. 
Thus, in making the adiabatic approximation we not only neglect higher order 
time derivatives, but we also neglect higher order crossed and time deriva- 
tives. We will see shortly, that, the solutions of equations I.52a and I.52b, tend 
to evolve slowly after some initial transients, but, they also tend to develop a 
large gradient in fj., in the vicinity of /x = 0. This large gradient invalidates the 
neglect of the term containing the crossed derivatives. In the next paper in this 
series an alternate to equation 1.55 will be developed in which the crossed de- 
rivative term is retained. We will see that this term plays an important role in 
approximating the solutions of equations I.52a and I.52b. 

The numerical solutions of the system of partial differential equations (1.44- 
1.46) were obtained by an explicit "marching" method on a uniform /j, -space mesh, 
but with a variable, self-adjusting time Step. We write the system 1. 44-1.46 in a 
condensed vector notation as 

P = K(Q) 1.56 
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where, 



1.57 


and, 




K(0) = 


g _ h + (1 - ^2) 


2^ 




\ - l/x| h + g 


1.58 


and then the shortest dynamical time scale in the system at the current time 


r" = ^ a/ is, 


r 

nni 


T “ min 


Qi 


i“1.3 < 

j=l,NPTS 


S"Ql 




1.59 


where the subscripts indicate components of the vector equation 1.57; the left 
superscript indicates the index associated with the time marching, and the right 
superscript indicates the index associated with the /i -space mesh, r" is the 
shortest time over which "substantial” changes occur in any one of the integra- 
tion variables, Qj. 

1 

Since we sought an accurate evolution for f(/^,T) = 3 N/B/j-, we augmented 
equation 1.59 by introducing 


T • r: 

' f min 


f j 


1.60 


j=l,NPTS |3fj 

,Br , 
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and then by defining the shortest evolutionary time scale, r^= minimum (t" 

The quadrature time step Ar", after n time steps, was then set to some conven- 
ient fraction of r"; for the quadratures discussed here, Ar"= 0.1 r". 

H/ Q E 

A threshold time step ^ was introduced so that node formation in any 

of the , or initialization (equation 1.49), could not cause Ar^ = 0 for any n. The 
threshold time step was adjusted by investigating the fraction of time when it 
superseded the quadrature time step, and by the effects its size had on external 
quality figures of the quadrature. Thus, the final form of our time step algorithm 
is given by ArJ = maximum (At^, ). 

It is clear from equations I.52a and I.52b, that at /i. = 0, C(0, r) is a true con- 
stant of the evolving solution which by equation 1.51 is given by, C(0,r) = N(0, 0). 
Thus, 

=N(0, 0) 1.61 

This condition is not a part of the explicit marching algorithm but is checked 
after each integration step as a measure of the fidelity of the overall finite 
difference scheme. 

In Figure 1 , we show an example of a plot of the left hand side of equation 
1.61 for the solution shown in Figure 5. Also shown is the reference constant, 
C(0,0). Note the departures of C(0,t) from the analytic constant after approxi- 
mately one Larmor period, coincident with the appearance of oscillations in 
the steep gradient in f(/ 2 , t) near /j. = 0. Beyond this point the finite difference 
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Figure 1. Evolution of C(0,r) (see equation 1.61) for the broad parallel beam injection 
shown in Figures. The failure of C(0 ,t) to remain constant coincides with the develop- 
ment of the erratic oscillations in Figure 5 in the vicinity of /i = 0. 




equations are unable to faithfully replicate the behavior of the continuum solu- 
tions, and the numerical solution is terminated. For all other numerical solu- 
tions presented in this paper, C(0, r) has been monitored and has been found to 
remain very nearly constant as in the first Larmor period of Figure 1. 

We have considered a variety of implementations of the Laplacian finite dif- 
ference operator with no perceptible change in the results reported here. The 
solutions which will be discussed shortly have generally been checked by doubling 
the mesh density in /x and finding no significant variations in the numerical results. 

A numerical quasi-linear adiabatic solution was obtained from equation 1.55 

by treating it as a standard diffusion equation with a spatially dependent diffusion 

coefficient. From equation 1.55 it is clear that N(0,r) = N(0, 0). The zero in the 

diffusion coefficient is equivalent to an impenetrable membrane at /X = 0. By 

placing a yU-space mesh point at ft = 0, we have prevented any "leakage" from one 

half-space to the other. The finite difference equation is implemented via the 

00 

explicit method of DuFort and Frankel which is unconditionally stable. 

All solutions reported graphically in this paper depict j+i _ 

- fJ) in a connect-a-dot fashion. No smoothing has been made in order to 
extract "f . The plots were made on a Calcomp plotter with 0.01" resolution 
with the distance between Larmor period tickmarks being 6.25 inches. 

In Figure 2, we show an explicit numerical solution of the quasi-linear adiabatic 
equation for an anisotropic initial condition corresponding to a narrow beam of 
particles injected with Gaussian probability about the mean magnetic field direction. 
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(IN GYROPERICX)S) 


Figure 3. Evolution of a narrow beam injected parallel to the mean mag- 
netic field in the quasi-linear diabatic approximation to the gyro-phase 
averaged probability distribution function. 


The three dimensional isometric presentation gives the solution through an 
elapsed time of two Larmor periods of the particle motion in the mean mag- 
netic field. In this solution, and all those to follow, = 0.09. Because of the 
zero in the diffusion coefficient of equation 1.55, probability does not propagate 
through yu. = 0; the step which invalidates the adiabatic approximation within the 
quasi-linear framework is clearly shown. The magnitude of this step grows in 
time until a uniform density is established through the forward pitch angle cone 
(as proven in Appendix B). 

The neglected terms in equation I.52a, which make the quasi-linear diabatic 
description different from the Markovian quasi-linear adiabatic description of 
equation 1.55 , ^ make important differences in the time evolution and asymptotic 
states of the probability distribution function. A numerical solution of equations 
I.52a and I.52b, for the same initial condition as that discussed in the previous 
paragraph, is presented in Figure 3. It is immediately clear that making the 
adiabatic approximation within the quasi-linear approximation, does serious in- 
justice to the evolution which it purports to approximate. The evolution of the 
probability distribution function proceeds smoothly through yU, = 0 , the site of the 
step formation of the adiabatic approximation. Furthermore, at the end of two 
Larmor periods, a substantial number of particles have passed through 90°. 

We notice that after two Larmor periods the diabatic solution (Figure 3) still 
contains a narrow strong- gradient transition region which is represented in the 
adiabatic approximation by the step. We expect this narrow transition region. 
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Figure 4. Evolution of a narrow beam injected across the mean magnetic 
field in the quasi-1 inear diabatic approximation to the gyre -phase averaged 
probability distribution function. 







and the formation of the corresponding step in the adiabatic approximation, for 
any initial condition for which N(0,0) 4 1/2. 

The evolution of the quasi-linear diabatic solution for cross-field injection, 
is illustrated in Figure 4. In this solution, the initial condition corresponds to 
a narrow beam of particles injected with Gaussian probability centered about the 
direction perpendicular to the mean magnetic field. Low frequency oscillations 
are clearly present at early times near /x = 0; the telegraph-like transients which 
reach the /x -space bounds at approximately 0,8 Larmor periods are rapidly 
damped. In the two Larmor periods shown, the probability distribution function 
attains a more nearly isotropic state than in the parallel injection case considered 
previously, and the narrow transition region which developed in that case is not ‘ 
apparent in this solution. As we will demonstrate, this absence of the transition 
region is due to the even (in /i) initial condition of this solution. 

One further quasi-linear diabatic solution is illustrated in Figure 5. In this 
solution, the initial condition corresponds injection of a broad beam of particles 
along the mean field direction. The ratio of probabilities in the forward to the 
backward directions in this beam is approximately nine to one with a very sharp 
gradient in probability distribution at = 0. This sharp gradient causes the 
immediate generation of a pair of oppositely propagating wave fronts which damp 
rapidly in the vicinities of /u. = ±1. For post transient times, this solution takes 
on a shape which is similar to that shown in Figure 3. In the next section of this 
paper, we will show that this similarity can be understood in terms of a 
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mode expansion of the probability distribution function. We have concluded that 
the high frequency oscillations which develop in this solution, in the vicinity of 
= 0, beyond approximately r = 1.0, are not a real feature of the solution. In 
the presence of the high frequency oscillations in Figure 5, the computed value 
of C(0,r) was found to oscillate considerably about its previously constant value. 
(See Figure 1 and the discussion surrounding equation 1.61.) We have further 
cheeked the invariance of all the diabatic solutions which we have presented here 
to changes in the/x -space grid spacing. The high frequency oscillations are not 
invariant to changes in grid spacing, but instead, manifest themselves in a variety 
of forms which depend on the grid spacing. In the particular solution presented in 
Figure 5, the ^-domain has been divided into two hundred intervals; every tick 
mark in the figure corresponds to two grid points. 

VI. A MODE EXPANSION OF THE PROBABILITY 

DISTRIBUTION FUNCTION 

In this section we introduce a ’’mode" expansion of the probability distribu- 
tion function which provides considerable insight into the post transient evolution 
of the quasi-linear diabatic solutions. Within the context of this mode expansion, 
the reduction of the quasi-linear solutions to the adiabatic limit is easily under- 
stood, and the reasons for the failure of the adiabatic limit become apparent. It 
is this point of view which has provided the primary motivation for the success- 
ful kinetic approximation which will be introduced in the next paper in this series. 
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The modes which are introduced in this section are obtained by Laplace 
transforming the quasi-linear diabatic equation, equation 1.20, and then finding 
the eigenfunctions generated by the Laplace transform of the integrodifferential 
operator on the right hand side of that equation. We call the eigenfunctions of 
the Laplace transformed operator "Laplace-modes.” This procedure is analogous 
to the standard modal analysis of plasma physics. 

The integral of the right side of equation 1.20 is a convolution integral under 
Laplace transformation. Thus, the Laplace transform of equation 1.20 is. 


p) 1.62 

9/i. 

where, 

Sfo. p) = (1 -/x2)KO, p) 1.63 


p?0x, p) - f (^, 0) = rp’ 

ajj. 


and K (Ai,p) is the Laplace transform of the kernel with p for the Laplace vari- 
able. Unless otherwise stated, only real, non-negative p will be considered here. 
The normal modes are introduced through 






’bfj. 


+ ^m(p) P) = 0 


1.64 


In Appendix A, we show that iK (/x ,p) is positive definite in p- for real, positive p. 
However, when p = 0, X {p,Q) takes on a zero at /x = 0. Thus, (/^,P) is non- 
negative in /X, with isolated zeroes at /x = ±1, and also at /x = 0 when p = 0. The 
situation v/hen p = 0 will be considered separately in a moment. When p 7 ^ 0, the 
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eigenfunctions which are the solutions of equation 1.64 are the solutions of the 
Sturm-Liouville problem. Thus, these eigenfunctions form a complete, orthog- 
onal set which can be used to form a mode expansion of f(/x, p). 

CO 

^0. p) = ^m(p) P) 

m *= 0 

By substituting this expression into equation 1.62 and using the orthogonality of 
the eigenfunctions, we find, 

J p)T(/x, 0) 

f.(P) = — 1.66 

EJ(P) [p + ’7^\(P)1 

where, 

K (p) ~ J (>> p) 1-67 

Equations 1.65 and 1.66 together form a representation of the exact solution of 
equation 1.20. 

One immediate solution of equation 1.64 follows from setting k = 0, and ^ = 
constant. This solution is p-independent. Furthermore, this eigenvalue has the 

Thus, we set Xq “ 

1.68 

The isotropic (a -independent, m - 0) part of the initial probability distribution 
function remains constant in time. 


2 3 

minimum value allowed; all other eigenvalues are positive, 
and =1. From equation 1.66, 

f»(p) = ^| 0). 
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In the limit, p = 0, a different situation presents itself. Because, (/j., 0) 
contains a zero at = 0 , the Sturm-Liouville method for the generation of the 
eigenfunctions cannot be applied as above. However, in the half domain, -1 </x < 0, 
where ^ (^,0) ^ 0, except on the boundaries, the Sturm-Liouville method can be 
applied. Because equation 1.64 is even in ii, this set of eigenfunctions is appro- 
priate for the positive half domain, 0 < /^ < 1, as well. For each half domain 
eigenfunction, two mutually orthogonal eigenfunctions in the full ^-domain can be 
constructed by taking even and odd combinations of the half domain eigenfunctions. 
With this choice of symmetry in [x, these eigenfunctions in the full -domain 
represent the limits, as p -* 0, of two of the (p ^ 0) eigenfunctions discussed 
above. Thus, in the limit, p = 0, the eigenfunctions which are the solutions of 
equation 1.64 become doubly degenerate; each eigenvalue corresponds to an even 
and odd pair of eigenfunctions. 

Notice from equation 1.62, if p is set to zero in ^ {/x,p), then it becomes 
the Laplace transform of equation 1.26. Thus, If (^,p) is replaced by J9 (fi,0) 
in equation 1.62, the adiabatic, or Markovian, approximation to the probability 
distribution function is generated. The eigenbasis for these adiabatic solutions 
can be generated from equation 1.64 by setting p = 0. Thus, the doubly degenerate 
eigenfunctions discussed above are the "adiabatic eigenfimctions"; they are 
reached by taking the limit, p = 0, in the "diabatic eigenfunctions" which are the 
solutions of equation 1.64 for positive p. Otherwise, the adiabatic eigenfunctions 
can be constructed from the even and odd combinations of the half domain 
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eigenfunctions as discussed above. An important property of these eigenfunctions, 
which is apparent from this construction, is that the odd adiabatic eigenfunctions 
are discontinuous at /x =0. 

As an example of the construction of the adiabatic eigenfunctions, consider 
the following simple case. Let domain eigenfunction 

which corresponds to this choice is (fi) = 1. The adiabatic eigenfunctions 
which can be constructed from this choice are, 

(ju, 0) = 1 1.69 

and 

(^) = xjj^ (^, 0) = 1 (/X < 0) 

= - 1 (/X > 0) 1.70 

We have already seen that 4>q is independent of p, so the choice given by 

equation 1.69 is obvious. The second choice, equation 1.70, follows from the con- 
siderations that ( /X, p) must be an odd function of [m with only one node in /x for 

any p; ^//(^^ (/x) is the only possible limit of 0i(fx,p). We will demonstrate that 
this choice for (/x ,0) is also consistent with the monotonic increase of the 
eigenvalues with increasing order at a fixed value' of p. Notice that 0) is 
discontinuous at /x = 0 as stated above. 

We have found that = 0, and from the Sturm-Liouville method, we can 
assure ourselves that all other adiabatic eigenvalues are greater than zero. We 
have concluded that X-j (p) - 0 as p 0, but, we can now further conclude that all 
higher order eigenvalues reduce to the positive adiabatic eigenvalues as p - 0. 
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With these results, we are able, in Appendix B, to explain the presence of the 
step in Figure 2 , as well as to prove the failure of the adiabatic approximation 
to propagate probability through i//^ = 0. 

In Figures 6 and 7, we present (/U-,P) and \p^{/j.,p) for several values of p 
ranging from zero to infinity. These curves were obtained using a Runge-Kutta 
integration routine and the "shot-gun method" on equation 1.64. The correspond- 
ing eigenvalues are listed in Table I. The adiabatic eigenfunctions were obtained 
by constructing the half-space eigenfunctions and then taking the appropriate com- 
binations of these as discussed above. The approach to the adiabatic limit is 
demonstrated in these figures for several small values of p. All of these numer- 
ical results were obtained for the exponential correlation function introduced in 
equation 1.43. In this case, 

3^0, p) = P ^ '.^1. - 1.71 


^ 1 



Since ^ (/i,p) is essentially independent of /x for very large p, the eigenfunctions 
are well approximated by the Legendre polynomials in this case, and the eigen- 
values are well approximated by, 

X (p) ^ 1.72 
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1 Table I 

1 

■j 

■f 

li 

1 The First and Second Eigenvalues for Various Values 


1 of the Laplace Variable, p 


i P . ; 

4 

CVJ 


GO 

0 

0 

10^ 

o.ooe 

0.005994 

10^ 

0.01992 

0.0596 

10^ 

0.1908 

0.5593 

1.0 

0.940 

2,66 

10-^ 

0.566 

2.64 

10-2 

0.275 

2.36 

lo'^ 

0.1705 

2.317 

lO-** 

0.2232 

2.322 

0 

0 

2.250 
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Thus, for non-zero m, we have, = 0, and the eigenfunctions plotted in fig- 

ures 6 and 7 for p = “ are just the Legendre polynomials with the appropriate 
normalization. 

We have plotted \ (p) and (p) in Figure 8 for the values of p given in 
Table I. The sol*d lines in that figure give various approximations to the p- 
dependence of the eigenvalues. Equation 1.72 has been used to obtain the solid 
line fits to the numerical results for large p. For small p, the solid line fits 
are given by 


\(P) 


\<A) 


-ln(p) 


1.73 


with Cj = 1.13 and C 2 = 0.571. We have been unable to construct an argument 
for equation 1.73 other than the observation that it works very well for small p. 

A visual comparison of Figures 3 and 5 with Figure 6, and of Figure 4 with 
Figure 7, has given us the impression that, on top of the contribution of (a^,P) = 1 
the long time behavior of the numerical solutions for is dominated by the 

lowest order eigenfunction which has a non-zero projection onto the initial condi- 
tions for the probability distribution function. In each case, we have the impres- 
sion that evolves, with increasing r, into distributions in A which have the 

appearance of the appropriate eigenfunctions with decreasing values of p. If this 
impression is correct, then a Idnetic approximation to the probability distribution 
function should be obtainable through any simplification of equation 1.62 which 
provides an accurate description of the low order eigenfunctions and eigenvalues 
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when p is small, but does not retain the full complexity of that equation when p 
is not small. In paper II of this series , we will construct a Idnetic approximation 
using this guiding principle; here, we present a crude argument which supports 
our visual impression. 

From Figures 6 and 7, we see that both</^j^ (f^.P) and are relatively 

independent of p for small p and jx not too close to zero. We will assume that 
the large time behavior of tip- ,r) is given by the low order eigenfunctions with 
small p. Therefore, we approximate equation 1.65 by, 

2 

P) = 0 ) 

m = 0 


with. 


fm(P) = 


(m. 0) fO. 0) 


£2(0) [p + (p)j 


1.75 


and then, 

2 

f Oi, '7") = 0). 1.76 

m“0 

First, consider the case m = 2. Under our working assumption, and for any 
f (/x,0) which is even in /x, the dominant time-dependent part of f (p.,r) for large 
r, should be given by (m> ^)* obtain an approximation to f 2 (r) by using 

equation 1.75 and also (p) = which can be justified from equation 1.73 for 

small p. After these approximations, we see that we expect £^(r) a 
i.e. , we 62 q)ect the non-isotropic part of the distribution function to decay away 
exponentially with a predicted decay rate. To test this prediction, we have studied 
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the numerical solution for f (/x , r) given in Figure 4. This solution has a Gaussian 
distribution in /i, centered at /x = 0, for an initial condition; in particular, the ini- 
tial condition is even in . The amplitude of f 2 (t") was determined by measuring 
the difference between the value of the numerical solution and the final isotropic 
level at /x ^±0.3, where, from Figure 7 we can see that (/u.»p) is nearly con- 
stant over a wide range in values of p. This amplitude is plotted in Figure 9 as 
a function of time, and compared to a strictly exponential decay with the predicted 
decay rate. (The amplitudes of both curves are arbitrary; they have simply been 
placed near each other for easy comparisons.) In spite of the level of this argu- 
ment, the agreement between the numerical solution and the exponential decay 
predicted here is quite good. The fact that the ?igreement in decay rate does not 
seem to persist to very large times is due to our inability to accurately measure 
i^{r) once its amplitude has become very small. 

If the initial probability distribution function is not even in , then we expect 
the long time distribution in /x to be dominated by the m = 1 term in equation 1.76. 
We do not predict an exponential decay in this case, however, since we cannot 
approximate (p) by a constant value for small p. Instead, we have (P)p 7 
Cj/(-ln(p)) and fi(p) ^ - ln(p) for small p. Erdelyi^'^ has constructed an Abelian 
theorem for Laplace transforms which states, if a function of time behaves like 

(lnr)“r^“^ r - ® 1.77 

then its Laplace transform b^- aves like, 

(-lnp)‘"p"^ p-0^ 1.78 
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AMPLITUDE (ARBITRARY UNITS) 


1.0 



Figure 9. The large time decay of the first and second eigenmodes 
compared to their Abelian predictions. 


5] 



where a = 1 is included, but k should be greater than zero. We will apply this 
theorem here for = 0 anyway, and find that it apparently makes some sense 
even in this limit. We have determined f^(r) in a manner similar to that out- 
lined above for f 2 (T), but we have used the numerical solution presented in Fig- 
ure 3, and have measured the amplitude of the non- isotropic part of f(Ai,T) at 
jj.- ±1. The amplitude of fj(T), so determined, is plotted in Figure 9 along with 
a fit to these results which is proportional to ln(r)/r. Once again, the fit between 
the time dependence of the measured amplitude and the predicted behavior is 
quite good. 

Although this argument is not conclusive, we feel that the evidence weighs 
ill the direction of our original supposition. The broad features of f seem 
to be dominated for large r by the low order eigenfunctions. Since fo('^) is a 
constant in time, and f ^(r) decays in time so slowly compared to the exponential 
decay of the higher order modes, we expect that a truncation of the Laplace-mode 
expansion beyond m = 1 would give a reasonable description of the average be- 
havior of f(/x,r) for large r. Of course, if f (/x, 0) happened to be an even func- 
tion of jj., then the m - 2 term in the mode expansion should also be included. 

VII. DISCUSSION 

One of the major purposes of this series of papers is to determine whether 
or not the quasi-linear approximation to the probability distribution function can 
be considered a leading approximation in a systematic expansion procedure. We 
have concluded in this paper that the presence of a strong mean magnetic field 
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can preclude the application of the adiabatic approximation to the quasi-linear 
probability distribution function. Thus, the Markovian limit, in which the propa- 
gation of a charged particle in magnetostatic turbulence is governed by a Fokker- 
Planck equation in velocity space, cannot be reached. As a corollary to the above 
conclusion, we further conclude that the quasi-linear adiabatic approximation to 
the distribution function cannot be used to judge the success or failure of the pro- 
posed quasi-linear expansion scheme. We will show in paper IE of this series 
that the quasi-linear approximation does contain an intrinsic non-unifonaity in/x . 
This defect is not related to the failure of the adiabatic approximation. 

We have proceeded with a study of the properties of the quasi-linear solutions 
with no further approximations in this paper. Numerical solutions of the quasi- 
linear diabatic equation have been obtained for an axisymmetric slab model of the 
plasma turbulence with an exponential two-point correlation function. Striking 
new wave phenomena have been discovered in the diabatic solutions. In these 
solutions, propagation through /z = 0 has been found, in contrast to the adiabatic 
solutions. On the other hand, this propagation through /x = 0 has been found to be 
very slow. For the case of a beam injection along the mean magnetic field direc- 
tion, the relaxation to isotropy has been found to proceed like ln(r)/(r) in contrast 
to the familiar exponential relaxation associated with diffusive behavior. 

A Laplace-mode expansion of the probability distribution function has been 
introduced. Because we have dealt in this paper with the quasi-linear diabatic 
equation which is integrodifferentlal in time, rather than differential with constant 
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coefficients, many of the familiar properties of the analogous modal analysis 
have not been readily apparent, but have been shown to still exist. In particular, 
it does seem that the long time evolution of the probability distribution function 
is dominated by the lowest order Laplace-mode which has a non-zero projection 
on the initial distribution function. It also seems that the long time evolution is 
determined by the low order Laplace-modes from the vicinity of the Laplace 
space origin. We have shown that ^ the origin (the adiabatic limit) the Laplace- 
modes become doubly degenerate and discontinuous functions of /z , thus leading 
to the failure of the adiabatic approximation. In the next paper of this series, we 
construct a new kinetic approximation to the quasi-linear solutions by consider- 
ing the vicinity of the Laplace space origin, rather than just the origin, thereby 
removing the degeneracies and discontinuities in the Laplace-modes. 
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APPENDIX A 


BASIC PROPERTIES OF THE LAPLACE TRANSFORMED 
KERNEL FUNCTION 


The quantity, K (yu-,p) is defined through, 


p) = dr R(jj-, r) cos r 


When p = 0, K(/x,0) = Mq (/x), where M„(/x) is the zero'th moment of the kernel 
whose properties have been discussed in the main text of this paper following 
equation 1.39. In particular, we have seen in the text that Mq (/u) is positive defi- 
nite except at /X = 0, where it is zero. In the following we consider the case p > 0. 

We introduce the Fourier transform of the correlation function (see equation 


1. 10) through. 


R(r) dk R(k) 


By introducing this expression into equation A.l, and then inverting the order of 


integration, we find, 


+ (k/i - 1)2 p2 + (k/x + 1)’ 
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If R(r) were a correlation in space or time which was obtained through space or 


time averaging in the following sense: 


R(r) = Urn 

L-OD 


-f 

21 - 1 , 


dz B'(z) B' (z + r) 


A.4 


then, by Cramer's law, R{k) would be non-negative for real In the case be- 
ing considered in this paper, in which the correlation function is defined through 
ensemble averaging, we assume that R(k) is non-negative in order to preserve 
the properties of the more realistic correlation functions. Then, from equation 
A.3, we see that when p > 0, the integrand which appears in that equation is non- 
negative for any fj- and k. Thus, we conclude thati^ (M.p) is positive definite for 
p > 0 and for any unless R(k) is zero for all k; i.e., unless the random field is 
set to zero everywhere. 
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APPENDIX B 


FAILURE OF THE ADIABATIC APPROXIMATION TO 
PROPAGATE PARTICLES THROUGH = 0 

The adiabatic eigenfunctions which were introduced in the main text in the 
discussion preceding equation 1.69, make up a complete set of orthogonal eigen- 
functions which can be used to construct the adiabatic approximation to the prob- 
ability distribution function. Using the method which lead to equations 1.65 and 
1.66, we find, 


in which, 


where , 


CO 

f O, p) = 2^ f,„(p) (P-) 

m“ 0 

f "f (/^. 0) 

‘'-I 


f (p) = 


[p + 


■r 

•^-1 




B.l 


B.2 


B.3 


The minimum eigenvalue which is available to us, from this construction, is 
\(A) = 0, and the corresponding eigenfunctions (doubly degenerate) are = 1, 

and ~ 1 fj- < 0 and 'Ai^^ = -1 for 0 < fx. All other eigenvalues are posi- 
tive, and part of a discrete spectrum. 

Because of the simple dependence of these equations on the Laplace variable, 
p, the Laplace inversion can be done. We find 
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in which, 


where, 


Clearly, 


CO 

f (/Li, r) = E 


B.4 


m ® 0 


f fr) = f (0) e 




B.5 


E(0)=— f 0) 

m A 


B.6 


f 0^, r) f; fo (0) i//(A) (^) + f ^(0) i/;(A) 5_7 

from which, we generally expect the adiabatic approximation to the probability 
distribution function to approach a final state which contains a discontinuous step 
at /Li == 0. 

We have also seen, in the discussion leading to equation 1.69, that the adiabatic 
eigenfunctions can be constructed from the half-domain eigenfunctions which form 
a complete, orthogonal set in either of the half- domains, -1 < /ll ^ 0, or 0 ^ M ^ 1. 
Consider the half-domain, -1 < /Li < 0. In this half-domain, both members of each 
degenerate pair of adiabatic eigenfunctions are identical, and equal to one of the 
half-domain eigenfunctions. Thus, each member of each pair of adiabatic eigen- 
functions is orthogonal, in the half-domain, to all other eigenfunctions which are 
members of other pairs. Therefore, from equation B.4, 
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Thus, the total probability in each of the half-domains is conserved in the adiabatic 
approximation; propagation through = 0 is impossible. It is important to note 
that the conclusions stated in this appendix depend only on the following properties 
of the power spectrum (see equations 1.38 through 1.41); P(w) must be non-zero 
for 1 < CO < ® , and 
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